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Abstract 

We rely on a variational approach to derive a set of equations governing a trapped 
self-interacting Bose gas at finite temperature. In this work, we analyze the static 
situation both at zero and finite temperature in the Thomas-Fermi limit. We derive 
simple analytic expressions for the condensate properties at finite temperature. The 
noncondensate and anomalous density profiles are also analyzed in terms of the 
condensate fraction. The results are quite encouraging owing to the simplicity of 
the formalism. 
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1 Introduction 



In a remarkable series of experiments on Rubidium and Sodium vapors|Tl [2], the Bose- 
Einstein condensation was first observed. Altliough having been predicted theoretically 
a long time ago for noninteracting boson systems[3], the experimental challenge was to 
demonstrate that a real gas can indeed be "bose condensed". Since then, a great effort 
was devoted by researchers all around the world in order to understand and predict the 
condensate properties. The main tools, beside the Monte-Carlo calculations[l], were the 
Bogoliubov[5], the Popov|6], the Beliaev[71 E] and the Hartree-Fock-Bogoliubov[9[ [TOl ITT] 
approximations. These approximations all adopt some simplifying assumptions about 
the various quantities involved in the problem, such as the order parameter or the 
condensate density Uc = the non-condensed density or thermal cloud n and the 

anomalous density m. A major well-known drawback of these methods is that they cannot 
be easily extended to situations where their main assumptions fail. In a previous paper [12]. 
we rely on a different approach, based on the time- dependent variational principle of 
Balian and Veneroni[13j , which allows one to overcome some of those restrictions. We 
obtained a set of three coupled dynamical equations, which we called "Time-Dependent 
Hartree-Fock-Bogoliubov" (TDHFB) equations, governing the evolution of $, n and rh. 
They were shown to generalize in a consistent way the Gross-Pitaevskii equation [14] . 

The present paper is devoted to analyze the implications of our TDHFB equations. 
Since it is important to apprehend first what happens in the static situation before going 
further into the analysis of the excitation spectrum or to the full dynamical case, we 
will focus on the static solution both at zero and finite temperature in the Thomas- Fermi 
(TF) limit. The interest is evident, since there remain many unanswered questions such as 
the general dependence of the density profiles on the temperature and on the interaction 
strength and the effect of the interactions on the critical temperature. More particularly, 
recent experiments are raising challenging questions about the precise determination of the 
thermal cloud and its backeffects on the condensate [15]. Indeed, due to the difficulties 
inherent to these experiments, there is no clear image on the way the condensed and 
non condensed phases mix up. Hence, in these preliminary calculations, we intend to 
provide some simple answers. We do not pretend of course to reproduce exactly the 
experimental data or the full Monte-Carlo calculations, but we would like to show that 
the simplifications that we are actually using (Mean field + Thomas- Fermi) are controlable 
and retain also the most important qualitative features without destroying the underlying 
physics. This provides a simple enough tool which can be considered as a starting point 
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for a more elaborate treatment. 

The paper is organized as follows. In section 2, we recall the main steps that have 
been used in [12] to derive the TDHFB equations. Then, we present the static solutions 
and discuss their properties at zero temperature. At finite temperature, the equations are 
much more involved and require a careful analysis. In the TF limit, we present a simple 
method which allows for a self-consistent determination of the various density profiles as 
well as some other static properties of the condensate such as the chemical potential and 
the condensate radius. Indeed, the TF approximation obviously provides simple enough 
analytical expressions since it neglects the kinetic terms thus yielding algebraic equations 
instead of partial differential equations. This is the main advantage of our method which 
yields the most important qualitative features without having to handle highly non-linear 
differential equations. 

In section 3, we present the results of our calculations. We plot first the condensate 
radius and the central density as functions of the condensate fraction and note in particular 
the compression effect of the condensate due to the thermal cloud. Moreover, we discuss 
the TF profile obtained for the condensate density even at low condensate fraction. The 
noncondensate density profile is also plotted for a wide range of condensate fraction and 
shows a good qualitative agreement with recent experiments. Finally, the anomalous 
density, although not yet measured experimentally, is shown to behave in a quite intuitive 
way. 

Some concluding remarks are given at the end of the paper. 



2 The TDHFB Equations and Their Static Solutions 

The general TDHFB equations were derived in ref . [12] for a grand canonical Hamiltonian 
of trapped bosons with quartic self- interactions (with coupling constant g and mass m): 

a(r) + - / a+(r)a+(r)a(r)a(r). (2.1) 
2 Jv 

The quantity V"ext(r) is the trapping potential and /i is the chemical potential. These 
equations read: 

^^•^ = (-|5^A + Kxt -yu + ^nc + 2^n) $ + c/m<l>*, 

ihfi = g (m*^^ - m$*2) , (2.2) 
ihm = g{2h + 1/V)<!>^ + 4 (-|^A + Kxt - /i + 2c/n + f (2n + 1/V)) m. 



H 



2m 



A + V;xt(r) -/i 
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where we have introduced the volume V of the gas in order to ensure the correct dimen- 
sions. In Eqs. fl2.2p . $ is the order parameter, ric the condensate density {ric = n the 
non-condensed density (or thermal cloud) and m is the anomalous density. The quantity 
n = ric + n is the total density. 

The TDHFB equations with a general Hamiltonian H were derived in [16]. The 
properties discussed here and in [12] were established there in their general forms. These 
equations were obtained using the Balian-Veneroni variational principle [13], with a gaus- 
sian trial density operator (that is, an exponential operator of a quadratic form) in the 
creation and annihilation operators. The result was a set of coupled evolution equations 
for the expectation values (a), (a+a) — (a+)(a) and (aa) — When one identifies 

these quantities respectively with the order parameter $, the non-condensed density n 
and the anomalous density m, and when one restricts H to the class ( 12. ip . the equations 
(1221) follow. 

The TDHFB equations couple in a consistent and closed way the three densities. 
They should in principle yield the general time, space and temperature dependence of 
the various densities. Furthermore, they obviously constitute a natural extension of the 
Gross-Pitaevskii equation[ll]. They are not only energy and number conserving, but 
also satisfy the Hugenholtz-Pines theorem (see below) which leads to a gapless excitation 
spectrum in the uniform limit. Moreover, the two last equations in (12.21) are not totally 
independent since n and m are related by the "unitarity" relation [12]: 

I = {l + 2Vnf -{2V\m\)\ (2.3) 

where the Heisenberg parameter I (which is always > 1) is a measure of the temperature, 
the lower limit being the zero temperature case. For instance, for a thermal distribution 
at equilibrium, / writes as / = coth^ {huJo/2kBT), where uq is the average frequency of 
the trapping field [12] . We therefore see that upon replacing n by its expression given 
in (12. 3p . the temperature appears explicitly in the equations. 

It is to be mentioned that the TDHFB equations have also been derived by several 
authors using different variational formulations |17l ITS] . In the first reference, the authors 
have obtained a set of equations very similar to ours. In fact, we can show that our 
equations can be deduced from theirs by taking the diagonal elements (r = r ) of the 
equations (Bl), (B2) and (B3) @ of ref.[n]. 

The static solutions, which are the object of our study in this work, are obtained 
by setting to zero the right hand sides of (12. 2p . At zero temperature, the standard TF 

^In fact, one can show that for a system of energy E, \fl =1 + 2 fsiE), where Jb is the Bose-Einstein 
distribution. 

■^There is however a factor 1/2 missing from (B3). 



4 



limit [19] amounts to neglecting the kinetic (or A) term in the Gross-Pitaevskii equation. 
This is particularly satisfied for strong interacting regimes or large atom numbers. At 
finite temperature and below the transition, since there are two phases (condensed and 
non condensed) which coexist, one has to provide a complementary recipe for what we 
shall call the finite temperature TF limit. First, neglecting the kinetic energy of the 
condensate remains a justifiable approximation since the atoms are slowed down in order 
to obtain condensation. On the other hand, m is believed to be an extremely small 
and slowly varying function whatever the temperature is (recall that it describes the 
correlations between the condensed and non-condensed phases). Hence, one may in a first 
approximation safely neglect Am. Heuristically, one may argue that, since the equations 
for ric and m contain almost comparable operators, and + g{nc + (1 + 2Vn)/AV), 
where /iq is the self-consistent mean field hamiltonian = Vext{r) — fJ' + gnc + "ign, the 
TF condition >> T {T being the kinetic operator), if fullfilled for ric shouldf also be 
satisfied for m. For this approximation to be consistent, and fh should vary on the 
same characteristic length, which is indeed the case as we will show later. 

Before proceeding further, it is important to notice at this point that a kinetic-like 
term of the thermal cloud does not appear explicitely in the equations but is rather hidden 
in the third equation of (12.21) . Indeed, the kinetic term of the thermal cloud is related to 
the second derivative of the anomalous density. Differentiating ( 12.31) yields a relation of 
the form: 

An ~ (V|m|)^ - (V|n|)^ + |m|A|m|, (2.4) 

which shows in particular that neglecting Am does not necessarily mean neglecting An. 
That is precisely the recipe that we shall adopt below. 

With this finite temperature prescription, the static equations corresponding to 
(12. 2p now write 

(Kxt(r) - fi + gn^ + 2gn) $ + gm^* = 0, 

m*$2-m$*2 = o, (2.5) 
(Kxt(r) -fx + 2gn+ ^(2Vn + 1)) m + ^{2Vn + l)<^^ = 0, 

These equations are naturally gapless and satisfy the Hugenholtz-Pines theorem [9]. In- 
deed, owing to the second equation in (12.51) . one can easily show that at zero momentum, 
the relation fi = g{n + n — \Th\) is clearly satisfied without adding further assumptions, 
as is usually performed [9]. 

In order to solve these equations, we may distinguish two rather different situations. 
The first one is for T = 0. When all the atoms are condensed, n = m = 0, and Uc equals 
the total density n of the gas. Omitting the trivial solution with Uc = 0, one may take into 
account just the first equation in (12.51) . since we consider a gas without a quantum cloud. 



5 



Indeed, within the present set of equations, it is an approximation (although justifiable) to 
ignore the quantum depletion at T = 0. The last two equations in (12. 5p become therefore 
meaningless, and we are left with a simple expression for the condensate density 

nc{r) = -e(r) = - (/i - Kxt(r)) . (2.6) 
9 

Upon defining the oscillator length oq = {h/mtuoY^'^ and the s-wave scattering length a = 
mg/Airh'^, we obtain for a spherical trapping potential Kxt('^) = ^rnujQr'^, the condensate 
radius R and the reduced chemical potential vq = ii/ \Tiujq for a gas of N bosons as 

1/5 

I , (2.7) 

2/5 



R 








ao 


V ao 


Z/Q = 







(2.8) 



The preceding expressions show that the spreading of the condensate depends essentially 
on the balance between the self-interactions and the trapping potential. These results 
have also been obtained by many other authors, see e.g. [HI [HI fT2] . 

When < T < Tbec; we have of course ric 7^ and ri 7^ 0. Let us introduce the 
parametrization 2Vh + 1 = \/7coshcr, 2l^|m| = -\/7sinho", which automatically endows 
the relation (12.31) . Then, from the third equation in (12.51) . one obtains a simple equation 
for X = e": 

3X^ - + 1 + 4Vne .2 _ 3W = 0, (2.9) 



from which one extracts n and |m| as functions of Uc- Next, one uses these expressions in 
the first equation (12.51) to get the condensate density 



^c(r) = -^(r) - - Vl-l]. (2.10) 




What is remarkable is that the sole acceptable solution of equation (12.91) is a bounded 
function oi t] = Vhc/Vi. It is represented on figure 1. 

Due to this behavior, one can easily show that the quantity which appears 

in (12.1 op is almost independent of ric and becomes rapidly close to unity. Indeed, since 
equation (I2.10p may also be rewritten as rj = — (X + 3X~^)/4, its solution provides 

the typical linear behavior shown in figure 2. 

Hence, one may safely approximate (I2.10p by 

n,(r)~-e(r)-:i(V7-l). (2.11) 

In fact, one can check that the relative error between the two expressions ( I2.10p and 
(12.111) is less than 1%. Finally, since vT does not depend on space, the result (12. lip 
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Fig.l: Solution of Eq.(2.9) vs. the dimensionless condensate density. Fig.2: Plot of the solution of eq.(2.10) showing a typical linear behavior. 



shows that the finite temperature correction to the Thomas-Fermi profile (12.61) is simply 
a space-independent (but temperature dependent) shift. This shift may be absorbed in a 
redefinition of the chemical potential which now writes 

/i = Kxt(i?) + -^(V7-l), (2.12) 

where R is the condensate radius. The condensate density finally writes in the suggestive 
form 

n.(r) = ^ext(i?)-Kxt(r) ^ ^2.13) 
9 

which is formally the zero temperature TF profile. It is then easy to show that the 
condensate radius takes also a simple form 

R / a \ 

^15nA] , (2.14) 



but now, it is the number of condensed atoms which is involved and not the total 
number of atoms. The same conclusion may be drawn for the chemical potential (I2.12p . 
Hence, our finite temperature prescription for the TF approximation provides natural 
extensions of the zero temperature expressions, since the Thomas-Fermi parameter is 
now Nca/ao instead of Na/aQ. 

In order to apprehend better these results, let us compute the remaining unkown 
quantities, such as the non condensed and the anomalous densities. To this end, and in 
order to obtain tractable expressions, we find it more convenient to use the simple fit 

X ^ (2.15) 

r] + 2/3 

(instead of the full analytical solution of equation (12.91) ) which reproduces correctly the 
solution X plotted in figure 1 with a residual error less than 0.1%. Upon rewriting 
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equation fl2.13p in the form t] = rio{l — x^), with an obvious definition of t]q, we obtain 
the non condensate density 



nix) 



■ 2V [ 2 \ rio{l - X 
and the anomalous density 



+ 2/3 v^r/o(l-x2) + 2/3, 



(2.16) 



(2.17) 



_l_^/J f V3rio{l-x^)+2/3 _ r7o(l-x^) + 2/3 \ 
2F 2 1, r/o(l - x2) + 2/3 V3r7o(l - x^) + 2/3/ ' 
as functions of the reduced radial distance x = r/R. In figures 3 and 4, we show typical 
density profiles (in reduced units) for tiq = 1. 




2Vmr 




x=r/R 

Fig.3: Noncondensate density vs. the radial distance. 



x=r/R 

Flg.4: Anomalous density vs. the radial distance. 



It is interesting to notice that, within our TF approximation, the figure 3 is quali- 
tatively consistent with the profile of the thermal cloud as depicted by Gerbier et a/|15j. 
It also predicts that the thermal cloud does not vanish at the boundaries of the conden- 
sate which is also compatible with the experimental results since it is widely known that 
condensed atoms are surrounded by the thermal cloud. What is less known (even exper- 
imentally) is the anomalous density (figure 4). This quantity behaves quite differently 
from the thermal cloud and is clearly dominated by a TF-like shape. In contrast with the 
thermal cloud, the anomalous density vanishes at the boundaries which is also plausible 
since the condensate vanishes there. Furthermore, we observe that ric, \rh\ and n vary on 
the same length scale (R) which justifies a posteriori our previous assumption. 

In order to obtain more quantitative results, one must determine Nc by using the 
normalization condition. We get easily the relation 

1 + 2N = 2N^ + Vlg{s), (2.18) 

where 



2 , N 3 , , 1 1 Is , 1 

q(s) = —p= + (v3 — l)s<l V s + larctanh , H — J —= + larctanh , 

W3 ' \ 2 2VV3 



(2.19) 
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with s = Ay/l /l5Nc. But since the function g{s) satisfies 1 < g{s) < 2/\/3, the equation 
(12.181) is approximately solved to yield, to a very good accuracy, the simple result 

N^^N- ^-^ — . (2.20) 

All the unknown quantities may now be determined in terms of and vT alone. The 
corresponding results will be discussed in the next section. 

3 Results and Discussions 

First of all, the condensate radius (I2.14p may be written as 

R = i?TF (^) , (3.1) 

where i?TF is the zero temperature result given by equation (12. 7p . Figure 5 represents the 
condensate radius (in units of -Rtf) as a function of the condensate fraction and we notice 
in particular the compression of the condensate when reducing Nc/N (that is increasing 
the temperature). This effect is by now a well established experimental result [15] and is 
attributed to the thermal cloud. The same effect of compression is observed on figure 6 
for the central condensate density nc(r = 0) but it is more pronounced due to the power 
law of 2/5 (see 12.131) instead of 1/5 for the condensate radius. To be more precise, let 
us choose generic values for the number of atoms and the interaction strength (A^ = 10^ 
and a/aQ = 0.5 10^'^) and plot the various densities (in units of the oscillator volume aj]) 
versus the radial distance (in units of Rtf = 3.758ao) for a condensate fraction ranging 
from 5% up to 60%. 




Fig. 5: Condensate radius vs. the condensate fraction. Fig. 6: Central density vs. the condensate fraction. 



The figure 7 shows typical Thomas-Fermi profiles for the condensate density, even 
for low condensate fraction. This is of course what one may expect on general grounds 
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in the TF regime. Moreover, the effect of compression of the condensate is also clearly 
visible here. 




The noncondensate density (from which we have substracted a constant n{R) for 
clarity) is plotted on figure 8 with the same units as before. The information which is 

NJN=5% 




1 1 1 1 1 1 1 1 1 1 

0,0 0.2 0,4 0.6 0,8 1.0 

Fig. 8: Noncondensate density vs. the radial distance for 
varying condensate fraction for #=10' and a/ao=0.5 10 '. 



added here with respect to figure 3 is the temperature dependence which appears via 
the condensate fraction. As noted earlier, the experimental result [T5] is well reproduced 
qualitatively for the whole range of condensate fractioijfl. In particular, we notice that 
when increasing the condensate fraction, the thermal cloud tends to spread and flatten. 

Although, in order to compare exactly with the experimental results, one must include an overall 
scale factor due to the finite ballistic expansion time. 
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On the other hand, the thermal cloud takes on a (small but) finite value for r > R. Even 
if this behavior is less intuitive, it is not very surprising since we do know that neglecting 
the second order derivatives amounts to making a cut of the densities at the boundaries. 
It is indeed a limitation of the TF approximation at the boundaries |20]. The tail should 
be reproduced when one reinjects the second derivatives in the equations. 




Fig. 9: Anomalous density vs.the radial distance for varying 
condensate fraction for A^=10-^ and a/uQ^O.S 10'-'. 



Finally, we plot on figure 9 the anomalous density. To our knoweledge, this quantity 
was never measured experimentally, and it is interesting that our calculations predict a 
very simple and yet intuitive behavior as well which remains of course to be confirmed. 



4 Concluding Remarks 

We present in this paper a finite temperature analysis of the static TDHFB equations 
(derived in a previous paper) in the Thomas- Fermi limit for a gas of bosons in a harmonic 
trap. 

At zero temperature, we obtain familiar expressions for the chemical potential and 
the condensate radius. The standard Thomas-Fermi profile for the condensate density is 
also recovered. 

At finite temperatures and below the transition, since there are two phases, one 
should provide a prescription for the TF limit. We propose such a recipe (maybe the 
simplest) which consists in neglecting the second order derivatives of the condensate den- 
sity and the anomalous density. The underlying idea is that, although the anomalous 
density is necessary for the coherence of the equations, it is believed to be a very small 
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and hopefully a very smooth quantity. We therefore obtain analytical expressions for the 
condensate density, the condensate radius, the chemical potential and the condensate frac- 
tion as functions of the temperature. Our expressions appear as natural extensions of the 
zero temperature TF limit, since the relevant parameter which controls the approximation 
becomes Nctt/ao instead of Na/aQ. 

Most importantly, we derive quite simple expressions for the noncondensate density 
and for the anomalous density, which we plot as functions of the condensate fraction and 
draw many conclusions. First of all, the condensate profile is almost of the TF shape of 
which the spatial extension and the heights are controlled by Nc/N. Furthermore, the 
compression of the condensate by the thermal cloud with increasing temperature is clearly 
visible. On the other hand, the noncondensate density profile is qualitatively consistent 
with the condensate fraction dependence observed in recent experiments. In particular, 
the thermal cloud tends to spread and flatten with increasing temperature. The calculated 
anomalous density, although not yet observed experimentally, shows also a very intuitive 
behavior; it is maximal at the center of the trap and zero at the boundaries. The tendency 
to spreading and flattening with increading temperature is also observed here. 

At the borders of the trap (where the condensate density vanishes) and for a given 
temperature, the non condensate density takes on a finite value which is a quite abrupt 
behavior. Athough this meets the fact that the thermal cloud is actually surrounding 
the condensate, it is to by no means conclusive. But this is also a shortocoming of the 
TF approximation as a whole since it is known to break down at the boundaries of the 
condensate. Indeed, reinjecting the second derivatives of the densities will certainly entail 
a more physical behavior. 

We are grateful to P. Schuck and Y. Castin for fruitful discussions and a careful 
reading of the manuscript. 

One of us (M. B.) is thankful to the members of the Groupe de Physique Theorique, 
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